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CONVERSION  FACTORS,  INCH-POUND  TO  METRIC  (SI) 
UNITS  OF  MEASUREMENT 


Inch-pound  units  of  measurement 

used  in  this  report 

can  be  converted  to 

metric  (SI)  units  as  follows: 

Multiply 

By 

To  Obtain 

cubic  feet 
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cubic  meters 

feet 
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meters 
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meters  per  second 

parts  per  thousand  (ppt)- 

0.02831685 

grams  per  liter- 

cubic  feet 

cubic  meters 

square  feet 

0.09290304 

square  meters 

square  feet  per  second 

0.09290304 

square  meters  per 

second 
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DEVELOPMENT  OF  A  NUMERICAL  SOLUTION  TO  THE  TRANSPORT  EQUATION 

Report  3:  TEST  RESULTS 

PART  I:  INTRODUCTION 

This  report  constitutes  the  third  report  in  a  three-report  series 
and  presents  the  development  of  a  test  program  and  the  numerical  experi¬ 
mentation  of  several  alternate  finite  difference  solution  techniques  for 
the  constituent  transport  equation.  The  following  schemes  as  previously 
developed  [1,2]  were  considered  initially  on  uniformly  spaced  grids. 

1.  Forward  Time  Upwind  Space  (FTUS) 

2.  Forward  Time  Centered  Space  (FTCS) 

3.  Spread  Time  Derivative  Upwind  Space  (STUS) 

4.  Spread  Time  Derivative  Centered  Space  (STCS) 

5.  Flux-Corrected  Transport  (FCT) 

A  7  x  8  uniform  grid  was  initially  employed  to  verify  the  coding  and  to 
evaluate  the  effectiveness  of  the  alternate  techniques.  Experimental  re¬ 
sults  are  directly  compared  with  known  analytical  solutions.  Analytical 
solution  techniques  are  initially  presented  followed  by  the  test  results. 
A  60  x  60  uniform  grid  was  subsequently  employed  to  test  boundary  condi¬ 
tion  effects  and  spatial  resolution  accuracy  relations. 

Based  upon  the  uniform  grid  test  results,  the  FCT  scheme  was  in¬ 
corporated  into  the  U.  S.  Army  Engineer  Waterways  Experiment  Station  Im¬ 
plicit  Flooding  Model  (WIFM) .  The  FCT  scheme  was  tested  for  a  sharp 
front  problem  over  the  variable  spaced  global  grid  employing  6783  com¬ 
putational  cells  representing  Mississippi  Sound. 
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PART  II:  ANALYTICAL  SOLUTION  METHODS 


The  advection,  diffusion,  and  the  combined  problem  are  considered 
here  for  the  constant  coefficient  case  in  turn. 

1.  Advection  Problem 

Consider  the  following  problems 


3c  3c  3c 

3t  "  -u  3x  "  V  3y 

(1.1) 

c(x,y,o)  =  <Kx,y) 


where 

c  =  constituent  concentration 

u  =  velocity  component  in  the  x  coordinate  direction 
v  =  velocity  component  in  the  y  coordinate  direction 
x,y  =  Cartesian  coordinates 
t  =  time 

Dc  3c  .  3c  ,  3c  .v  , 

Assume  u,v  are  constant.  Since  ^-“-j^+u  —  +  v  —  »  U*l)  “ay  be 
written 


Dc 

Dt 


0 


(1.2) 


From  this  relation,  we  observe  that  the  property  of  the  fluid  is 
unchanged  with  respect  to  a  coordinate  system  attached  to  the  fluid 
particles.  We  thus  expect  the  following  solution,  c(x,y,t)  ■ 
i|»(x  -  ut,y  -  vt)  .  Since 


3c  ,  3[4>(x  -  ut,y  -  vt)]  ,  }  3[^(x  -  ut,y  -  vt)].  {  * 

3t  3x  1  '  3y  V  ’ 


(1.3) 


our  suppositions  are  indeed  substantiated.  If  we  consider 


<p(x,y) 


x,y  e  (-R/2.R/2) 
x,y  t  (-R/2.R/2) 


as  our  initial  distribution,  then 


c(x,y,t) 


(-R/2  +  ut ,R/2  +  ut) 
(-R/2  +  vt ,R/2  +  vt) 
(-R/2  +  ut,R/2  +  ut) 
(-R/2  +  vt,R/2  +  vt) 


' l>(x  -  ut,y  -  vt)  (1.4) 


The  distribution  is  translated  in  space  with  no  distortion. 

2.  Diffusion  Problem 

Consider  the  following  problem,  using  the  same  notation  as  in  (1.1) 


c(x,y,o)  *  i//(x,y) 


(2.1) 


We  follow  an  approach  developed  by  Polzhiy  [3]  for  the  one-dimensional 
equation.  Assume  a  solution  of  the  following  form: 

I(x,y,t)  =  c(x,y,t) 


Gr  _  2  2  \/?_22f 

I  e  a  a  cos  ax  da  J  f  /  e  a  a 


(2.2) 


cos  ay  da) 


Since 


(2.2)  satisfies  (2.1)  above.  Next  we  compute 


OO  /  00 

•  f  2  2,  /  r  2  2^ 

31  I  -a  at.  jlf-aat 

—  *  -a  I  e  sin  ax  da  I  I  e 

— OO  \  —00 


cos  ay  da)  (2.3a; 


OO  J  CO 

aT  r  2  2.  /  r  2  2 

31  I  -a  at.  ,  /  I  - a  a  t 

=  -a  I  e  sin  ay  dal  I  e 

—00  \  —00 


cos  ax  da)  (2.3b] 


We  next  integrate  the  first  factor  in  (2.2)  by  parts;  i.e.,  I udv  = 

-)  o  J 


f  -a2a2t 

uv  -  Ivdu  .  Letting  u  =  e  and  dv  *  cos  ax  da  ,  we  obtain 


OO 

/ 


2  2  .  2  2 
-a  a  t  ,  -1  .  -a  a  t 

e  cos  ax  da  =  x  sin  axe 


(2.4: 


+  2a2tx 


i  r  2  2 

1  I  -a  a  t 
I  ae 


sin  ax  da 


Noting  the  first  term  on  the  right  hand  side  of  (2.4)  is  zero  we  obtain 
the  following  alternate  form  for  (2.2)  using  (2.3a): 


i(x.y.t) 


2  • 
-2a  tx 


(OO 

•/ 


2  2 

-a  a  t 


sin  ax  da 


■) 


(2.5: 


x 


cos  ay 


-2a2tx 


-1  31 
3x 


By  integrating  the  second  factor  in  (2.2)  by  parts  and  using  (2.3b)  we 
also  obtain  a  similar  relationship  with  x  and  y  interchanged. 


I(x,y,t)  =  -2a2ty 


sin  ay  da 


(2.6) 


^  J  e  a  3  C  cos  ax  da^  =  -2a2ty-1  ~ 


Equations  (2.5)  and  (2.6)  suggest  the  following  relation  for  I(x,y,t)  : 


f  -  I  £  -  I  -  Kx.y.O  -  Ce-(xV)/<^2t)  (2.„ 

3x  2aZt  3y  2a2t 


Consider 


I(o,o,t)  = 


Q  e-°2a2t  do^  ^  j  e-“2“2t 


dal  =  C 


(2.8) 


It  may  be  shown  that  each  factor  is  equal  to  -  . 

7T  a/t 

Thus  C  -  — 2~  and  we  finally  obtain: 
a  t 


Kx.y.t)  -  -§-  e'<*V>/(4a20 

azt 


(2.9) 


The  following  6  function  is  also  a  solution 


6(PfP',t-t')  -  I(x  -  x’,y  -  y',t  -  t') 


1 _  -r2/[4a2(t-t’; 


4a  ir(t-t') 


where 


r  -  J (x  -  x’)2  +  (y  -  y')2  ,  t  >  t' 


The  above  set  of  functions  with  parameter  a  -  t  -  t'  ,  have  the  three 


special  properties. 

00  oo 


f  f  4<P,P',C  -  t')  dx  dy  -  f - 1 -  e-(x-x') 

LL  L  - 1’) 

OO 

\  2a/n(t  -  t*) 


2/[4a2(t-t')] 


dx 


- (y-y ' ) 2/ [4a2 ( t-t ' ) ]  . 

e  dy  -  1  (2.10a) 


00 

f  j |«(P,P*,t  - 


t  * )  I  dx  dy  =  1. 


(2.10b) 


Urn  f  f 

t+t*  /  1  |6(P,P\t  -  t')|  dx  dy  =  0 


(2.10c) 


—CO  —  OO 


The  prime  on  the  integral  in  (2.10c)  denotes  integration  over  the  en¬ 
tire  x-y  space  except  for  an  arbitrarily  small  neighborhood  about  P'  . 
Thus,  after  concepts  of  Polzhiy  [3],  the  limit  element  of  the  sequence 
of  functions  generated  when  a+0  ,  in  the  sense  of  weak  convergence,  is 
a  6-function  with  respect  to  M*  ,  the  set  of  all  continuous  bounded 
functions  of  x  and  y  .  Thus  noting  6(P,P',t  -  t')  =  6(P',P,t  -  t'): 


lim 

t-*-t' 

+ 


1/ 


f(x') 


2aA(t  -  t') 


e-(x-x')2/[4a2(t-t’)]dx,( 


(2.11) 


11 


g(y') 


2a A (t  -  t') 


r<y-y')2/(4«2(t-f)ldyJ  ,  £(x)g(y) 


where 


ip(x,y)  -  f (x)g(y) 


11 


.T’.vr.v 


We  therefore  obtain  a  solution  of  the  form,  where  t'  ■  0 

c(x, 


00 

:,y,t)  =  /  f(x')  — i— 

J  2a/nt 

m  — OO 


e-(x'-x)2/(4a2t)  dx, 


[/  s‘y,)  ^ 


e-(y’-y)2/(4a2t)  dy, 


From  (2.11)  we  note  (2.12)  satisfies  the  initial  condition. 
Consider  the  following  problem 


f  (x)g(y)  -  <Kx,y) 


x,y  e  (-R/2,R/2)j 


o  x,y  t  (-R/2.R/2), 


and 


3c  /s2c  32c\ 

“  vv) 


Evaluating  (2.12)  we  obtain 


c(x,y,t) 


R/2  /Z~  2 

f  o  e“(x'-x)  /(^Kt)  dx, 

*  2V^irt 


L-R/2 


R/2 

/ 

L-R/2 


/r 


2i/Krrt 


o_  „-(y'-y)  / (4Kt)  dy, 


(2.12) 


(2.13) 


(2.14) 


Consider  (2.14)  in  the  following  manner: 


r  -  (x  -  x,)/2v'ScF  dr  -  -dxV2*€t 
s  -  (y  -  ds  -  -dy'/2*1ct 


then 


c(x,y,t) 


Noting 


(x+R/2)/2v£t  2 

e  r  dr 


-f  I 

*  (x-R/2)/2i4cF 


ft-  (y+R/2)/2»1ct 

—  / 

"  (y-R/2)/2(^t 


-s 

e  ds 


(2.15) 


/•  2  ?■ 
erf  (x)  =  —  /  e  r  dr  =  —  /  e  r 

AJ  A  J 

O  —X 


dr  , 


the  above  result  may  be  further  simplified  to  the  following  form: 

c(x,y,t)  -  jerf[(x  +  R/2)/2*1ct]  +  erf[(R/2  -  x)/2/Et] j 

x  |erf[(y  +  R/2)/2*£t]  +  erf[(R/2  -  y)/2^t]| 

3.  Advection  and  Diffusion  Problem 
Consider  the  following  problem 


Dc  3c  3c  3c  J 32c  .  32c\ 
Dt  "  3t  U  "  V  3y  / 


(2.16) 


(3.1) 


c(x,y,o)  =  <Kx,y)  =  f(x)g(y) 


where  u  and  v  are  constant.  Define 


G(x,y,t) 


{l£<x'fe 


-(x-x')Z/(4Kt) 


1  dx' 


{/8<y'fe 


-(y-y’)2/(4Kt) 


1  dy' 


(3.2) 


We  assume  a  solution  c(x,y,t)  -  G(x  -  ut,y  -  vt,t)  .  It  may  be  shown 
by  computing  3c/ 3t  ,  3c/ 3x  ,  3c/ 3y  ,  32c/3x2  ,  32c/3y2  ,  that  the 
solution  assumed  satisfies  (3.1).  Furthermore,  from  previous  work  (3.2) 
satisfies  the  initial  condition  and  are  assumed  solution  is  the  solution. 

We  consider  <Kx*y)  to  be  as  given  previously,  namely,  our 
familiar  block  distribution  problem.  Employing  previous  results,  the 
solution  assumes  the  following  form. 


c(x,y,t) 


jerf[(x  -  ut  +  R/2)/2/Kt)  +  erf[(ut  -  x  +  R/2)/2i1ct]  j 
*  {erf [(y  -  vt  +  R/2)/2^t]  +  erf[(vt  -  y  +  R/2)/2^Ct]  1  (3.3) 
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PART  Ills  SMALL  GRID  (7x8)  TEST  RESULTS 


Several  tests  were  performed  employing  this  grid  both  on  the 
Kirtland  Air  Force  Base  (CRAY-1)  and  the  Boeing  Computer  Service 
(CRAY-1)  computers.  The  problems  considered  and  numerical  results  are 
presented  in  turn.  In  characterizing  the  test  cases,  the  following  di¬ 
mensionless  parameters  are  employed: 


UT 

K  T 

dD  -  ~ 

Pe 

■v 

Ax 

x  ,  2 

Ax 

x 

"«* 

VT 

_  K  T 

D  x 

Pe 

y 

V 

Ax 

y =  7T 

Ay 

'tv 

where 

t  »  Simulation  time  step  length  (seconds) 

Ax  ■  x  space  step  (feet*) 

Ay  -  y  space  step  (feet) 

*  Dispersion  coefficient  in  the  x  direction 

Ky  ■  Dispersion  coefficient  in  the  y  direction 
V  ■  Particle  Courant  number  in  the  x  direction 

*  Particle  Courant  number  in  the  y  direction 

^  ■  Diffusion  number  in  the  x  direction 

■  Diffusion  number  in  the  y  direction 
Pex  ■  Peclet  number  in  the  x  direction 

Pe  ■  Peclet  number  in  the  y  direction 

y  7 

*  A  table  for  converting  the  inch-pound  units  of  measurement  found  in 
this  report  to  metric  (SI)  units  can  be  found  on  page  5. 
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1.  Block  Translation  Problem 

The  pure  translation  of  a  square  block  of  concentration  1  ppt 
initially  occupying  grid  cell  (3t6)  as  shown  in  Figure  1  was  considered. 
Particle  Courant  numbers  of  0.2  in  both  coordinate  directions  were  main¬ 
tained  for  ten  time  steps  of  100  seconds  duration. 

The  problem  represented  an  extremely  severe  test  for  the  solution 
schemes  considered,  since  a  discontinuity  was  introduced  over  only  1 
grid  cell  in  both  directions.  The  analytical  solution  consists  of  a 
unit  concentration  in  cell  (5,4)  after  100  seconds. 

Results  for  the  FTUS  scheme  are  presented  in  Figure  2.  The  scheme 
exhibits  considerable  frontal  smearing  but  is  nonosclllatory.  The  FTCS 

_3 

scheme  results  in  Figure  3  contain  122  versus  74  ppt  *  10  of  the 
FTUS  scheme  in  cell  (5,4)  at  the  end  of  the  simulation.  However,  the 
FTCS  develops  severe  oscillations  behind  the  front.  Results  for  the 
STUS  scheme  are  given  in  Figure  4  and  are  similar  to  those  of  the  FTUS 
scheme.  Results  for  the  STCS  scheme  of  Figure  5  are  similar  to  the  FTCS 
scheme  results.  The  FCT  scheme  results  in  Figure  6  are  nonosclllatory 
and  also  exhibit  a  certain  amount  of  frontal  smearing. 

All  schemes  considered  exhibited  numerical  dispersion  and  subse¬ 
quent  loss  of  mass  through  the  grid  boundary.  The  schemes  were  ranked 
in  Table  1  on  their  ability  to  maintain  material  on  the  grid.  The  FCT 
scheme  is  the  most  accurate  scheme  based  upon  this  criterion.  In 
Table  II  the  schemes  are  ranked  according  to  their  ability  to  place 
the  most  material  in  the  correct  cell  (5,4).  The  FCT  scheme  ranked 
first  in  maintaining  material  on  the  grid  and  was  also  the  top  ranked 
nonosclllatory  scheme  in  terms  of  mass  placement. 
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Figure  1.  Pure  translation  of  a  block  (small  grid  case) 

Time  Step  No.  10  Time  -  1000.00  Concentration  Field  (xio3) 
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Figure  2.  Pure  translation  of  a  block  (small  grid  case)  FTUS 
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Figure  3.  Pure  translation  of  a  block  (small  grid  case)  FTCS 
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Time  Step  No.  10  Time  =  1000.00  Concentration  Field  (xlQ3) 
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Figure  4.  Pure  translation  of  a  block  (small  grid  case)  STUS 
Time  Step  No.  10  Time  •»  1000,00  Concentration  Field  (xlO3) 
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Figure  5.  Pure  translation  of  a  block  (small  grid  case)  STCS 
Time  Step  No.  10  Time  ■  1000.00  Concentration  Field  (xlO3) 
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Figure  6.  Pure  translation  of  a  block  (small  grid  case)  FCT 
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Table  1.  Block  Translation  Scheme  Rankings  (Mass  Retention) 


Scheme 

7  3 

Material  Remaining  on  Grid  (ppt  x  10  ft  ) 

1.  FCT 

0.920873 

2.  STUS 

0.837683 

3.  FTUS 

0.812556 

4.  FTCS 

0.308442 

5.  STCS 

0.024073 

Table  II. 

Block  Translation  Scheme  Rankings  (Mass  Placement) 

Scheme 

Material  in  Cell  (5.4)  (ppt) 

1.  STCS* 

0.398 

2.  FTCS* 

0.122 

3.  FCT 

0.091 

4.  FTUS 

0.074 

5.  STUS 

0.072 

*  Oscillatory  solution. 

2.  Stationary  Block  Problem 

The  stationary  block  problem  specified  Initially  a  1-ppt  (%) 
concentration  in  grid  cell  (3,6)  as  shown  in  Figure  7.  Mo  diffusion  and 
no  advection  were  considered  and  the  schemes  were  tested  on  their  ability 
to  maintain  the  unit  concentration  in  the  grid  cell  for  10  time  steps 
of  100  seconds  each. 
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Figure  7. 

Stationary  block  problem  (small  grid  case) 

The  STUS  and  STCS  schemes  are  equivalent  for  this  problem.  The  re¬ 
sults  of  the  STCS  scheme  are  shown  in  Figure  8.  Note  the  tendency  to¬ 
ward  oscillation.  The  FTUS,  FTCS,  and  FCT  schemes  are  equivalent  for 
this  problem.  The  results  of  the  FCT  scheme  are  shown  in  Figure  9  and 
show  no  tendency  toward  oscillation. 

3.  Block  Diffusion  Problem 

The  straight  diffusion  of  a  square  block  of  1-ppt  concentration 
initially  occupying  grid  cell  (3,6)  as  shown  in  Figure  10  was  considered. 
Diffusion  numbers  of  0.0001  in  both  coordinate  directions  were  main¬ 
tained  for  ten  time-steps  of  100  seconds  each.  Analytical  solution  re¬ 
sults  Indicated  a  negligible  amount  of  diffusion  of  material  from  grid 
cell  (3,6). 

The  results  of  STCS  and  STUS  schemes  are  equivalent  for  this  prob¬ 
lem.  The  STCS  results  shown  in  Figure  11  exhibit  oscillation,  with 
7  3 

0.009  *  10  ppt-ft  of  material  being  diffused  from  cell  (3,6).  The 
results  of  the  FTUS,  FTCS,  and  FCT  schemes  are  equivalent  for  this  prob¬ 
lem.  Results  shown  in  Figure  12  for  the  FCT  scheme  exhibit  no  oscilla- 

7  3 

tion  and  diffuse  only  0.004  *  10  ppt-ft  of  material  from  cell  (3,6). 


m 


.  *  •  .  •*.  «* . 
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Time  Step  No.  10  Time  ■  1000.00  Concentration  Field  (xio3) 


N 


M 

1 

2 

3 

4 

5 

6 

7 

1 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

2 

0. 

-0. 

0. 

-0. 

0. 

-0. 

0. 

3 

0. 

0. 

-0. 

0. 

-0. 

0. 

0. 

4 

0. 

-0. 

0. 

-0. 

0. 

-0. 

0. 

5 

0. 

0. 

-0. 

0. 

-0. 

0. 

0. 

6 

0. 

-0. 

1000. 

-0. 

0. 

-0. 

0. 

7 

0. 

-0. 

-0. 

0. 

-0. 

0. 

0. 

8 

0. 

0. 

0. 

0. 

0. 

0. 

"7 

0. 

o 

(STUS,  STCS)  /Mass  Grid  Total  1.00000  x  107  ppt-ft3 


Figure  8.  Stationary  block  problem  (small  grid  case)  STUS  and  STCS 
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FTCS,  AND  FCT)  Mass  Grid  Total  1.00000  x  10  ppt-ft 
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Figure  10.  Diffusion  of  a  block  (small  grid  case) 


3 

Time  Step  No.  10  Time  =  1000.00  Concentration  Field  (*10  ) 
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(STUS,  STCS)  Mass  Grid  Total  1.00088  *  10  ppt-ft 
Figure  11.  Diffusion  of  a  block  (small  grid  case)  STUS  and  STCS 
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Time  Step  No.  10  Time  =  1000.00  Concentration  Field  (xlO  ) 


1 

2 

3 

4 

5 

6 

7 

0. 

0. 

0. 

0. 

0. 

0 

0. 

0. 

0. 

0. 

0. 

0 

0. 

0. 

0. 

0. 

0. 

0. 

0 

0. 

0. 

0. 

0. 

0. 

0. 

0 

0. 

0. 

1. 

0. 

0. 

0. 

0 

0. 

1. 

996. 

1. 

0. 

0. 

0 

0. 

1. 

0. 

0. 

0. 

0 

0. 

0. 

0. 

0. 

0. 

0. 

0 

(FTUS,  FTCS,  FCT)  Mass  Grid  Total  0.999999  x  10  ppt-ft 

Figure  12.  Diffusion  of  a  block  (small  grid  case) 

FTUS,  FTCS,  and  FCT 

4.  Front  Translation  Problem 

Constituent  concentration  in  cell  (1,4)  was  initially  set  at 

1  ppt  and  maintained  at  that  level  for  10  time-steps  of  100-sec  dura¬ 
tion.  The  velocity  was  directed  in  the  horizontal  at  a  strength  of 

2  fps.  No  diffusion  was  considered.  The  problem  setup  is  shown  in 
Figure  13. 

7  3 

Over  the  10  time-step  run  2  x  10  ppt-ft  of  mass  is  injected. 

7  3 

Since  each  cell  represents  10  ft  and  the  velocity  magnitude  is  2  fps, 
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Figure  13.  Front  translation  problem  (small  grid  case) 

cells  (1,4),  (2,4),  and  (3,4)  should  all  be  at  1  ppt  strength  at  the 

end  of  the  simulation. 

The  results  for 

the  FTUS  scheme  are  shown  in  Figure  14.  The 

Time  Step  No.  10  Time  -  1000.00  Concentration  Field  (xio3) 

Mass  Grid  Total  2.97691  x  10  ppt-ft 
Figure  14.  FTUS  front  translation  results 

7  3 

frontal  smearing  extends  to  the  boundary  of  the  2  x  10  ppt-ft  mass 

7  3  7  3 

Injected  and  the  1  x  10  ppt-ft  initial  mass,  2.97691  *  10  ppt-ft 

remains  in  the  systems,  the  remaining  mass  is  lost  through  the  boundary. 

The  results  of  the  FCT  scheme  are  shown  in  Figure  15.  The  frontal 

7  3 

smearing  is  less  severe  for  these  schemes.  Of  the  2  x  10  ppt-ft 
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Figure  15.  FCT  front  translation  results 


7  3  7 

mass  injected  and  the  1  x  10  ppt-ft  initial  mass,  2.99515  x  10  ppt- 

3 

ft  remains  in  the  system,  the  remaining  mass  is  lost  through  the 
boundary . 

5.  Rotation  Problems 

Consider  the  uniform  rotation  of  strength  w  (radian/sec)  about 
cell  (no,mQ)  as  shown  in  Figure  16.  Define  the  following  relations. 


x'  =  (n  -  nQ)Ax 


y'  -  (m  -  mQ)Ay 


where 

m,  n  =  indices  for  cell  (n,m) 

Ax,  Ay  =  grid  spacings  in  the  x  and  y  directions,  respectively 

m  ,  n  =  indices  for  rotation  center 
o  o 

For  an  arbitrary  cell  (n,m),  the  velocity  v  is  computed  as  follows. 

__  ^  ^  A  A  A  A  A 

v  1  u  x  r  =  uk  k  (x’i+y'j)  *  -wy’i  +  wx'j 


v 


* 

«,  * 

Figure  16.  Rotation  problem  notation 
Block  rotation  problem 

A  block  of  material  initially  occupying  cell  (3,3)  at  strength 
1  ppt  is  rotated  through  90°  over  10  time  steps  of  100  second  duration. 
At  the  end  of  this  simulation,  the  block  should  occupy  cell  (3,5)  at  a 
strength  of  1  ppt.  The  problem  setup  is  shown  in  Figure  17. 

The  results  of  the  FTUS  are  presented  in  Figure  18.  The  distri¬ 


bution  is  considerably  smeared,  with  mass  being  lost  through  the 
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Figure  17.  Rotation  block  problem  (small  grid  case) 

Time  Step  No.  10  Time  =  1000.00  Concentration  Field  (xlO"*) 
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Figure  18.  FTUS  block  rotation 
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boundaries.  Of  the  original  1  x  10  ft  -ppt  of  mass,  0.80215  x  io  ft  - 
ppt  remains  on  the  computational  grid. 

The  FCT  scheme  results  are  presented  in  Figure  19.  The  distribu¬ 
tion  is  again  smeared,  with  mass  being  lost  through  the  boundaries.  Of 

.  7  3  7  3 

the  original  1  x  10  ft  -percent  of  mass,  0.878427  x  io  ft  -percent  re¬ 
mains  on  the  computational  grid.  Thus,  the  FCT  scheme  is  superior  to 
the  FTUS  scheme  in  this  regard. 
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Figure  19.  FCT  block  rotation 


Front  rotation 

Constituent  concentration  in  cell  (1,4)  was  initially  set  at 
1  ppt  and  maintained  at  that  level  for  10-time  steps  of  100-sec  dura¬ 
tion.  A  circular  velocity  field  at  the  center  of  cell  (1,7)  was  sped' 
fled  such  that  the  front  would  rotate  and  occupy  cell  (4,7)  at  the  end 
of  the  simulation.  The  problem  setup  is  shown  in  Figure  20. 


Figure  20.  Front  rotation  problem  (small  grid  case) 


Over  the  10  time-step  run  4.7124  *  10  ft  -ppt  of  mass  is  in- 

7  3 

jected.  Since  each  cell  represents  10  ft  and  considering  boundary 
input  cell  (1,4),  the  total  grid  mass  should  equal  5.7124  x  107  ppt-ft 
The  results  for  the  FTUS  scheme  are  shown  in  Figure  21.  Frontal 


smearing  extends  to  the  boundary  resulting  in  a  grid  mass  total  of 
4.82762  x  io7  ppt-ft3. 


The  results  of  the  FCT  scheme  are  shown  in  Figure  22.  Frontal 


smearing  is  less  severe  in  this  scheme  with  a  mass  grid  total  of 
5.22026  x  107  ppt-ft3. 
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Figure  21.  FTUS  front  rotation 
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Figure  22.  FCT  front  rotation 


6.  Summary  of  Small  Grid  Test  Results 

By  employing  a  7  *  8  grid,  it  was  possible  to  verify  the  computa¬ 
tional  schemes.  Due  to  the  grid  size  boundary  effects  influenced  the 
computations  for  the  majority  of  the  problems  considered. 

The  FCT  method  proved  to  be  the  superior  technique  for  the  block 
translation  problem,  since  the  mass  lost  through  the  boundary  was  less 
for  this  scheme  than  for  all  other  schemes. 

For  the  block  diffusion  problem,  simulation  results  were  the  same 
for  FCT,  FTCS,  and  FTUS  schemes.  This  result  is  expected,  since  in  the 
absence  of  advection  the  FTCS  and  FTUS  schemes  are  equivalent.  The  STCS 
scheme  exhibit  minor  oscillations  in  the  solution  surface.  This  result 
is  somewhat  unexpected. 

In  both  the  block  and  front  rotation  problems,  the  FCT  scheme 
proved  superior  to  the  FTUS  scheme  in  terms  of  minimizing  mass  lost 
through  the  computational  boundaries. 


PART  IV:  LARGE  GRID  (60  *  60)  TEST  RESULTS 


A  60  x  60  computational  grid  was  constructed  to  perform  additional 
testing.  All  work  in  this  phase  was  conducted  on  the  Boeing  Computer 
Services  CRAY-IA  in  Renton,  Washington.  The  problems  considered  and 
scheme  results  are  presented  in  turn  below. 

1.  Single  Cell  Block  Translation 

The  same  problem  considered  in  the  small  grid  setting  was  setup 
as  shown  in  Figure  23.  The  initial  unit  concentration  in  cell  (30,30) 
was  to  be  advected  to  cell  (32,28)  over  ten  time  steps  of  100  second 
duration. 

The  results  for  the  FCT  scheme  are  shown  in  Figure  24.  We  note 
the  symmetry  about  the  diagonal  (the  direction  of  translation) .  All 
mass  is  preserved;  i.e.,  no  mass  is  lost  through  the  boundaries. 

2.  Single  Cell  Block  Diffusion 

The  diffusion  of  a  unit  concentration  initially  in  cell  (30,30) 
is  considered  as  shown  in  Figure  25  for  10  time  steps  of  100  second  dura¬ 
tion.  Diffusion  coefficients  are  increased  from  1  ft/sec  of  the  previ- 

2 

ous  small  grid  to  100  ft  /sec. 

The  results  of  the  FCT  scheme  are  shown  in  Figure  26.  Since 

-  Ky  ,  the  results  are  symmetrical  in  the  x  and  y  axes.  No  mass  is 

lost  through  the  boundaries.  The  computed  results  are  compared  with  the 

analytical  solution  results  in  Table  III.  The  computed  results  represent 

6  2 

values  over  an  entire  grid  cell  area  of  10  ft  ,  whereas  the  analytical 
results  correspond  to  the  point  value  at  the  center  of  each  cell.  As 
At, Ax, Ay  0  ,  the  computed  results  should  approach  the  analytical  value. 
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Figure  24.  FCT  block  translation  (large  grid  case)  (*10  ) 


Table  III.  Single  Cell  Block  Diffusion  Simulation 
Versus  Analytical  Solution  Comparison 


Location 

Analytical 

Result 

Computed 

Result 

(30,30) 

0.7768 

0.684 

(31,30) 

0.1018 

0.068 

(32,30) 

0.0003 

0.003 

(31,29) 

0.0172 

0.007 

3.  Multicell  Block  Translation 

For  this  problem  Ax  *  Ay  =  100  ft,  and  At  =  10  seconds  as  shown 
in  Figure  27.  100  cells  were  used  to  represent  the  same  distribution 

represented  by  one  cell  in  the  single  cell  problem.  In  order  to  reduce 
computational  costs,  50  time-steps  were  employed. 

The  FCT  scheme  results  at  the  end  of  50  time-steps  are  presented 
in  Figure  28.  Note  again  the  symmetry  about  the  streamline.  No  mass 
was  lost  through  the  boundary  and  no  oscillations  were  present  in  the 
solution.  All  material  in  the  multicell  block  translation  problem  is 
contained  within  an  area  equal  to  that  of  6  cells  of  the  single  block 
translation  simulation,  while  in  that  simulation  the  original  distribu¬ 
tion  was  spread  over  19  cells.  Thus  by  reducing  the  time-  and  space- 
steps  a  more  accurate  result  is  obtained. 

The  STCS  scheme  results  at  the  end  of  50  time-steps  are  presented 
in  Figure  29.  Oscillations  developed  behind  the  front  and  spread  to  the 
boundary,  essentially  destroying  the  solution.  Mass  appeared  to  be 
created  through  the  first  37  time-steps  and  destroyed  over  the  last 
13  time-steps. 
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Figure  28.  FCT  block  translation  (large  grid  case)  (xl(r) 


4.  Multicell  Front  Translation 


For  this  purpose  the  spatial  step  is  100  ft  and  the  time  step  is 
10  seconds  as  shown  in  Figure  30.  Constituent  concentration  in  cells 
(1,26) -(1,35)  was  set  initially  at  1  ppt  and  maintained  at  that  level 
for  50  time-steps.  The  velocity  was  directed  in  the  horizontal  at 
2  fps.  No  diffusion  was  considered. 

5  3 

Over  the  50-time-step  run  100  *  10  ppt-ft  of  mass  is  injected. 

5  3 

Since  each  cell  represents  10  ft  ,  the  sum  of  the  cell  concentrations 

3  3 

at  the  end  of  the  simulation  should  equal  110  ppt-ft  (100  ppt-ft  due 
to  injection  and  10  from  the  boundary  condition).  The  results  of  the 
FCT  scheme  shown  in  Figure  31  completely  preserve  mass.  Since  the  hori¬ 
zontal  velocity  is  2  fps,  at  the  end  of  the  simulation  all  concentra¬ 
tions  in  the  rectangular  area  (1 ,26)-(l,35) ,  (11,26)-(11,25)  should  be 
equal  to  1  ppt.  As  shown  in  Figure  31,  the  FCT  scheme  exhibits  some 
frontal  smearing. 

This  problem  demonstrated  the  ability  of  the  FCT  scheme  to  pre¬ 
serve  mass  for  injection  (source)  type  problems. 

5.  Summary  of  Large  Grid  Test  Results 

The  FCT  scheme  preserves  mass.  No  mass  was  lost  through  the  com¬ 
putational  boundaries  in  either  the  single  advection  or  diffusion  prob¬ 
lems.  The  multicell  advection  (block  translation)  problem  demonstrated 
that  by  reducing  the  time  and  space  steps,  the  FCT  scheme  will  produce 
a  more  accurate  result.  The  STCS  scheme  results  for  this  problem  ex¬ 
hibited  severe  oscillations  behind  the  front  eventually  destroying  the 
solution.  The  multicell  front  translation  problem  demonstrated  the 
ability  of  the  FCT  scheme  to  handle  injection-  (source-)  type  problems. 
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Bs8&d  upon  those  Additional  tests  in  conjunction  with  the  smell 
grid  test  results,  the  FCT  scheme  was  selected  for  incorporation  in 
the  Waterways  Experiment  Station  Implicit  Flooding  Model  (WIFM) . 
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Figure  30.  Multicell  front  translation  (large  grid  case) 
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PART  V:  VARIABLE  GRID  TESTING  ON  THE  MISSISSIPPI 
SOUND  GLOBAL  GRID 

The  FCT  algorithm  was  incorporated  as  a  group  of  subroutines  in 
the  WIFM  to  simulate  salinity.  Density  coupling  was  not  considered. 
Since  the  numerical  investigation  of  Mississippi  Sound  is  concerned 
only  with  fixed  boundary  problems,  the  flooding  routine  was  removed 
from  WIFM.  The  resulting  hydrodynamic-salinity  code  (WIFM-SAL)  enables 
the  treatment  of  the  general  fixed  boundary  nonlinear  problem  on  a 
variably  stretched  grid.  The  exponentially  stretched  grid  shown  in 
Figure  32  employing  6785  computational  cells  with  a  minimum  spatial 
resolution  of  approximately  4000  feet  has  been  developed  to  describe 
global  circulation  and  horizontal  salinity  variation  in  Mississippi 
Sound. 

In  order  to  describe  the  dispersion  mechanics  within  Mississippi 
Sound  the  following  relations  are  employed  in  WIFM-SAL. 

K  -  D^-^+  K’ 
x  C  x 

(5.1) 

Ky.D^M>  +  K. 

where 

2 

K  ,  K  =  dispersion  coefficients  (ft  /sec) 
x  y 

D  =  dimensionless  constant  (5.93-20.2)  in  the  direction  of 
flow  0.23  in  the  direction  perpendicular  to  flow 
C  =  Chezy  coefficient  (dimensions  of  /g) 
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g  =  gravity  (ft/sec  ) 

h  =  water  depth  digitized  over  the  global  grid  (ft) 
u,  v  =  velocity  components  in  the  x  and  y  directions, 
respectively  fps 

o 

K' ,  K '  =  additional  dispersion  effects  (ft  /sec) 
x  y 

As  in  the  testing  on  the  regularly  spaced  grids  the  following 

sharp  front  problem  was  considered.  The  hydrodynamic  problem  setup 

employed  initial  water  surface  elevations  set  to  zero  over  the  compu- 

2 

tational  domain,  flow  inputs  of  4  ft  /sec  per  unit  flow  width  at  cells 
(97,3)  and  (2,44),  and  a  ramp  function  of  1/6  ft  per  time-step. 

The  salinity  problem  setup  employed  initial  zero  concentration 
levels  over  the  entire  computational  region  except  in  a  block  of  ten 
cells  as  shown  below. 

n  t  (101,105) 
m  i  (54,55) 

n  e  (101,105) 
m  e  (54,55) 


Salinity  levels  were  maintained  at  zero  at  the  two  flow  input  locations. 

2 

In  the  dispersion  relations,  c  =  c  ■  10.0  and  K'  =  K'  =  0.0  ft  /sec. 

x  y  x  y 

A  5  time-step  simulation  was  performed  employing  a  time-step  length  of 
6  minutes  for  several  numerical  schemes. 

In  order  to  characterize  the  transport  characteristics  of  these 


simulations  the  following  dimensionless  numbers  may  be  considered. 


£ 


Cr 


u  . - /0  At 
n,mfl/2 


n,m 


Cr 


k  ] 

Vn+l/2 ,m| 

|  At 

n,m 


(li2)nM2 


(5.3) 


Pe 


k 


n,mfl/2 


n,m 


K 


xn,m-KL/2 


Pe 


k 

^n,m 


n+l/2,m 


(V2)„4a2 


K 


yn+l/2,m 


(5.4) 


where 


^  cell  (n,m)  particle  velocity  Courant  number  at 

*n,m 

time  k 


„  K 

Cr  =  x  -  a1 


Cr 


k 


y  =  y  -  ® 

n,m 


i2  cell  (n,m)  particle  velocity  Courant  number  at 
1. 


time  k 

_  k 

Pe  =  x  -  a, 

n,m 

cell 

(n,m) 

Peclet 

number  at 

time 

„  k 

Pey  =  y  -  «2 

n,m 

cell 

(n,m) 

Peclet 

number  at 

time 

k 

k 


At  =  time-step  length 
Aa^  =  space  increment 


A«2  =  a2 


u 


i2  space  increment 

(y^)m  =  cell  (n,m)  stretching  coefficient  in  the  direction 

(u2)n  =  cell  (n,m)  stretching  coefficient  in  the  ot2  direction 
1c 

n*4-i  /•>  =  velocity  component  in  cell  (n,m)  at  time  k  in  the  a, 
.rni-i/z  direction  1 

vn+l/2  m  =  velocity  component  in  cell  (n,m) 

*  direction 

Kxn  iiH-1/2  5  disPersion  coefficient  in  cell  (n,m)  at  time  in  the 
*  0^  direction 

k 


at  time  k  in  the  ct„ 

1 


JL, 

Kyn+l/2,m  B  Aspersion 
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V*  «*. 


In  general,  these  dimensionless  numbers  vary  in  time  as  well  as  in  space 
over  the  computational  grid  and  time  interval  of  concern.  Normal 
practice  is  to  replace  the  time  dependency  of  the  cell  velocity  compo¬ 
nents  by  their  maximum  values,  thus  removing  the  time  dependency.  In 
all  simulations,  the  following  relations  for  the  cell  particle  velocity 
Courant  numbers  hold. 


Crx  *  Cry  1  °-3 
n,m  yn,m 


(5.5) 


The  time  dependency  in  the  Peclet  numbers,  may  be  removed  by 


substituting  (5.1)  into  (5.4)  with  K'  *  K'  «  0.0  ft  /sec.  Thus  obtain 

x  y 

k  k 

for  Pe  (results  for  Pe  are  analogous) . 

Xn,m  ^n,m 


Vrt/zPl).*"! 


c  (p,)  Act, 
n,m  1  m  1 


D^h 


(5.6) 


!c  k 

Since  c  and  h  vary  extremely  slowly  with  time,  we  may 
n,m  n,m  J 

drop  the  k  superscript  and  obtain 


In  all  simulations,  the  following  relations  hold  for  the  cell  Peclet 
numbers  in  the  vicinity'  of  the  sharp  front: 

Pe  ,  Pe  >  100  (5.8) 

X  V  — 

n,m  7n,m 
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Results  at  the  end  of  the  simulation  for  the  FTUS  scheme  are 

shown  in  Figure  33.  All  concentrations  are  positive  and  the  initial 

14  3 

mass  of  0.132435  *  10  ppt-ft  equaled  the  final  mass  plus  the  diffu¬ 
sion  of  material  through  the  boundaries  to  within  the  precision  limits 
of  the  CRAY  I-S. 

Results  at  the  end  of  the  simulation  for  the  FTCS  scheme  are  shown 
in  Figure  34.  Since  the  cell  Peclet  number  limit  of  2  is  violated  in 
the  vicinity  of  the  front,  oscillations  develop  behind  the  movement  of 
the  front.  Mass  is  conserved  in  the  simulation,  at  the  expense  of  nega¬ 
tive  concentrations  and  cell  concentrations  greater  than  10.0. 

The  three  FCT  limiters  developed  previously  [2]  were  tested. 

The  original  Zalesak  limiter  results  are  shown  in  Figure  35.  Cell  con¬ 
centrations  greater  than  10.0  were  developed.  The  alternative  one 
(mixed  time  level)  limiter  results  are  shown  in  Figure  36.  No  cell 
concentrations  exceed  10.0.  However,  the  maximum  allowable  mass  to 
enter  a  cell,  Q+  ,  or  leave,  Q  may  now  be  negative  unlike  in  the 
original  limiter.  Some  small  negative  concentrations  are  also  devel¬ 
oped.  The  second  alternative  (previous  time  level)  limiter  results 
are  shown  in  Figure  37.  No  cell  concentrations  exceed  10.0,  Q+  and 
Q  are  nonnegative,  and  no  negative  concentrations  are  developed. 

Although  all  these  versions  of  the  limiter  preserve  mass,  the 
second  alternative  (previous  time  level)  limiter,  based  upon  the  non¬ 
negativity  of  Q+  and  Q  ,  appears  to  have  some  advantages  over  the 
other  two  limiters.  Additional  prototype  simulations  must  be  performed, 
however,  before  the  final  version  of  the  limiter  may  be  determined. 


Figure  37.  Alternative  two  (previous  time  level)  FCT 
limiter  global  grid  simulation  results  at.  5x  (x  10^) 
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